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Abstract 

Electron tomography is becoming an increasingly important tool in materials science for studying 
the three-dimensional morphologies and chemical compositions of nanostructures. The image quality 

'NOTICE: this is the authors version of a work that was accepted for publication in Ultramicroscopy. Changes resulting 
from the publishing process, such as peer review, editing, corrections, structural formatting, and other quality control mecha- 
nisms may not be reflected in this document. Changes may have been made to this work since it was submitted for publication. 
A definitive version was subsequently published in Ultramicroscopy (2013), DOI: 10.1016/j.ultramic.2013.01.002 
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obtained by many current algorithms is seriously affected by the problems of missing wedge artefacts 
and non-linear projection intensities due to diffraction effects. The former refers to the fact that data 
cannot be acquired over the full 180° tilt range; the latter implies that for some orientations, crystalline 
structures can show strong contrast changes. To overcome these problems we introduce and discuss 
several algorithms from the mathematical fields of geometric and discrete tomography. The algorithms 
incorporate geometric prior knowledge (mainly convexity and homogeneity), which also in principle 
considerably reduces the number of tilt angles required. Results are discussed for the reconstruction of 
an In As nanowire. 

1 Introduction 

Missing wedge artefacts and non-linear projection intensities due to diffraction effects are known to cause 
severe difficulties in electron tomography (ET) reconstructions obtained by standard methods. This has 
been reported, e.g., in [TH [TH1 EH ESI 133 SOI S3 [53] ■ Nevertheless, standard methods, such as filtered 
backprojection, algebraic reconstruction techniques, and simultaneous iterative reconstruction techniques [271 
[32] , are still widely used due to an apparent lack of alternatives [44] [15] . 

However, alternatives exist in the mathematical literature. Geometric tomography [17j . for instance, is 
concerned in part with the tomographic reconstruction of homogeneous (i.e., geometric) objects. Similarly, 
discrete tomography [281 129| usually deals with objects for which atomicity is a constraint or objects that 
exhibit a small number of attenuation coefficients. In many applications, certain prior knowledge about 
the shape of the structure of interest is available. For example, when reconstructing nanorods, nanowires 
or certain types of nanoparticles, one can typically assume that the structures are convex. (A subset K 
of points in the plane is convex if for any two points in K the line segment joining these two points lies 
completely within K; see Fig. [2]) In particular, in our experimental application of reconstructing an InAs 
nanowire from high-angle annular dark-field scanning transmission electron microscopy (HAADF STEM) 
data, it is known that the object is comprised of cross-sections that are mostly close to regular hexagons; 
see Q3J [52]. 

Here we demonstrate the use of geometric prior knowledge to overcome the problems of missing wedge and 
non-linear projection intensities due to diffraction effects by introducing four algorithms. For now, we use 
their abbreviated names; they are introduced in the next section. One of the algorithms (2n-GON) appears 
here for the first time and uses the strongest geometric prior knowledge available in our setup, namely that 
the slices contain nearly regular hexagons. Two algorithms, GKXR and MPW, are introduced here for 
the first time in the ET context and another algorithm (DART) is applied here for the first time to the 
reconstruction of a nanowire. As a fifth method, we discuss the BART algorithm, which was introduced in 
the 1970s and performs very well on homogeneous objects. BART has been implemented in the open-source 
software SNARK09 |12j and we provide commands and parameters that yield high quality reconstructions 
in our context. 

The idea of using geometric prior knowledge in ET appears already in [461 149] . In particular, |49] is, to 
the best of our knowledge, the only paper in ET that discusses geometric 2D slice-by-slice reconstruction 
methods, all of which are variants of the unfiltered backprojection (U-FBP) algorithm. Our work takes 
these investigations a step further and introduces alternative reconstruction methods, some of which are 
mathematically proven to converge towards the solution as the number of noisy measurements tends to 
infinity. We compare these methods (including U-FBP and the standard method SIRT as sixth and seventh 
algorithms, respectively) and show that these methods perform differently depending on the experimental 
setup and type of noise present in the data. 

The geometric reconstruction algorithms use prior knowledge to address the problem of the missing 
wedge and non-linear projection intensities. We show that reconstructions that are accurate (mean of the 
reconstructions is close to the true object) and precise (reconstructions have small variance) can be obtained 
by using data from considerably fewer tilt angles than in conventional tomography. This has practical 
consequences, because rapid data acquisition allows time resolved studies and imaging of beam-sensitive 
samples. 

As already mentioned, we investigate here only the use of geometric prior knowledge. Other types of 
prior knowledge might be related to the sparsity of the signal that represents the image gradient [24] [48] , or 
it might be assumed that the object is a realization of a random process with a given probability distribution 
[5] [TU]. For further applications of special-purpose reconstruction methods in a crystallographic context, see 

[Furjuns]. 
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2 Algorithms 



In this section we briefly describe the simultaneous iterative reconstruction technique (SIRT), the binary 
algebraic reconstruction technique (BART), the discrete algebraic reconstruction technique (DART), the 
Gardner- Kiderlen X-ray (GKXR) algorithm, unfiltered backprojection (U-FBP), the modified Prince-Willsky 
algorithm (MPW), and 2n-GON. Further details on SIRT, BART, DART, GKXR, and MPW can be found 
in [23l [26l [7j HH [20j 09]. Our results are based on Matlab implementations of these algorithms; BART is 
part of the open-source software SNARK09 . 

We consider only 2D versions of each algorithm, so that 3D reconstructions are obtained by 2D slice-by- 
slice reconstructions. Alternative reconstruction principles exist. For 3D and 2.5D approaches employing 
generalized Kaiser-Bessel window functions (blobs) for not necessarily homogeneous or convex objects, see 
|37M43j . We assume in this paper an acquisition geometry with a single tilt axis and a limited angular range. 

The algorithms require different input data. While SIRT, BART, DART, and GKXR take the projections 
as input, only the shadows are used in U-FBP, MPW, and 2n-GON. Note that following (non-mathematical) 
standard convention, projection refers to the measured intensity data (i.e., line integrals), while shadow 
denotes their support (i.e., the detector pixel locations that record non-zero intensities). It can therefore 
be expected that U-FBP, MPW, and 2n-GON are rather insensitive to intensity-affecting noise if the signal 
can still be distinguished from the background. Also note that the object's shadows represent binary data, 
which initially need to be extracted from the projection data. In this paper we achieve this by applying a 
suitable threshold and filter to the projection data; see Sections 13.21 and 14.21 For more advanced variants, 
such as edge enhancement, see [4"9"] . 

Another difference between the algorithms is that they use different types of prior knowledge. Most of the 
algorithms exploit the object's convexity or homogeneity. However, 2n-GON uses the additional assumption 
that the object is nearly a regular 2n-gon. A summary is given in Table [T] 
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Table 1: Overview of algorithms. 



2.1 Pixel-based reconstruction methods 

The first three methods discussed in this paper aim to reconstruct the individual pixel values of an image 
representing the object. 

2.1.1 Simultaneous iterative reconstruction technique (SIRT) 

In this paper we describe and employ an additive variant of SIRT [35J Section 7] . The additive SIRT algorithm 
(from now on referred to as the SIRT algorithm) is a standard technique for reconstructing grayscale images 
from tomographic data. When reconstructing homogeneous objects, an additional segmentation step is 
required that yields a binary image. We recall that SIRT is an iterative reconstruction algorithm that 
computes an approximate solution of the linear system Ax = b, where the vector x G ffi" contains the gray 
level (also referred to as pixel value) for each pixel, the vector b G R m contains the measured projection data, 
and the matrix A = ((%■) G ]g> m x™ represents the projection operation (i.e., computing the product Ax yields 
the projections corresponding to the image x). If no exact solution of this system exists, SIRT computes a 
solution for which the norm of the difference \ \Ax — b\ | (referred to as projection error) between the computed 
projection and the measured data is minimal with respect to a weighted L2 norm, i.e., a weighted least- 
squares solution (see [221 [25] for details). In contrast to other popular iterative algorithms, such as ART 
(algebraic reconstruction technique |27l Chapter 11]), the SIRT algorithm computes the projections for all 
angles in each iteration. Then the difference between these projections and the measured projection data is 
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computed. Subsequently, each image pixel value is updated by adding a weighted average of the projection 
difference for all lines that intersect this pixel. In other words, 

x (k+l) = x (k) _ XCA T D{Ax (k) _ 

where C = diag(l/ci, . . . , l/c„) £ R IlXTl with Cj = Y™=i a t] and D = diag(l/di, . . . ,l/d m ) 6 R mXm with 
di = Y^j=i a ij! the parameter A € K is the relaxation parameter of the algorithm (for our particular choice, 
see Sections 13.21 and 14.21) . 

2.1.2 Binary algebraic reconstruction technique (BART) 

The BART algorithm was introduced by Herman Along with other methods, BART is implemented in 
the open-source software SNARK09 [12] ■ The general idea is to enforce binary constraints on the solution x 
during the iterations of a chosen ART routine. In a second step, the solution is filtered to exclude isolated 
pixels. A complete code that can be read into SNARK09 (cither manually or as an input file) to obtain the 
BART reconstructions discussed in this paper is provided in Appendix A (Table [1} . 

2.1.3 Discrete algebraic reconstruction technique (DART) 

The DART algorithm, which has recently been proposed as a reconstruction algorithm for electron tomog- 
raphy [6l [7] is another algebraic method, in which a set of fixed pixels is updated in each iteration, reflecting 
the position of the boundary in the current reconstruction. The variant described here does not apply any 
subsequent filtering. For simplicity, we describe the DART algorithm for the specific purpose of reconstruct- 
ing a binary image. The general algorithm can deal with more than two gray levels. More on DART can be 
found in GSHEHiSI. 

Here, we focus on reconstructing a single object of homogeneous composition, resulting in a binary image 
reconstruction problem. The DART algorithm uses a continuous algebraic reconstruction method, such 
as SIRT, as a subroutine. From this point on, we will refer to this continuous method as the algebraic 
reconstruction method (ARM). After an initial gray level reconstruction has been computed using the ARM, 
this gray level reconstruction is segmented by global thresholding with threshold p/2, where the gray level 
p of the object is assumed to be prior knowledge. In practice, an appropriate value of p can often be 
obtained by first computing a SIRT reconstruction and then taking the average gray level over a region 
deeply in the interior of the object. One of the principal assumptions behind the DART algorithm is that 
errors in this segmentation are typically located near the boundary of the structure of interest. Indeed, 
when reconstructing a homogeneous object that is large with respect to the image pixel size, pixels that 
are deeply inside the interior of the object (e.g., a nanoparticle) are usually segmented correctly, while the 
segmentation of the boundary can be highly inaccurate, in particular when the ARM reconstruction suffers 
from missing wedge artefacts. The boundary can be computed from the initial segmentation as the set of 
pixels for which not all neighboring pixels belong to the same segmentation class. After the segmentation 
step, the set of pixels is separated into three subsets: the interior pixels I, the background pixels B and 
the boundary pixels F. Prior knowledge about the gray level p of the object and of the gray level of the 
background (here assumed to be 0) is now incorporated by solving the following constrained reconstruction 
problem, again using the ARM: 

solve Ax = b 
subject to 

Xi = 0, for i G B, 
x-i = p, for i 6 7. 

In this reconstruction problem, all interior pixels and background pixels are fixed to their respective gray 
levels and pixel values are only allowed to change for boundary pixels. This constraint strongly reduces the 
number of unknowns in the equation system, while the number of equations (i.e., the number of entries in 
b) remains unaltered. If the initial segmentation is of sufficient quality, the reconstruction of the boundary 
will significantly improve compared to the initial ARM reconstruction. 

The resulting reconstruction is again segmented, resulting in a new partition of the image into interior, 
background, and boundary pixels. As new pixels can be added to the boundary, pixels whose values were 
fixed in the previous step can now become boundary pixels and vice versa. The procedure of alternating 
segmentation and reconstruction steps is then iterated until a pre-defined convergence criterion is reached. 

A well-known limitation of DART is that the result can depend sensitively on the choice of gray levels 
and this choice may not be correct based on the first SIRT iterations. Sophisticated algorithms for choosing 
the correct gray levels can be found in [8j [50] . 
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2.2 Object-based reconstruction methods 



The second category of reconstruction methods that we consider in this paper is object based in the sense 
that the routines aim to determine a small number of parameters that completely describe the object (in 
our case, the vertices of polytopes). 

2.2.1 Gardner-Kiderlen X-ray (GKXR) 

The Gardner-Kiderlen X-ray algorithm [TH] is a recent development from the field of geometric tomography. 
It arose from theoretical work [22] in which it was shown that there are certain sets of four directions in 
2D such that the exact projections of a 2D convex object in these directions determine it uniquely among 
all 2D convex shapes. For example, directions specified by the four vectors (0, 1), (1,0), (2, 1), and (—1,2) 
constitute such a set [18] . The GKXR algorithm is based on the simple observation that given a sufficiently 
dense set of lines meeting a convex set K, the convex hull of all the points at which the lines intersect the 
boundary of K will form a convex polygon that approximates K well. The algorithm attempts to find this 
polygon for the set of projection measurement lines. 

Fig. Q] shows a schematic diagram of the basis of the algorithm. The unknown object is the oval K, 
assumed to lie inside the circle. For clarity, only a single projection, taken in the direction u, is considered in 
Fig. [2 although in practice projections in four different directions are used. For each projection direction u, 
detector pixels are located at the equally spaced points ti,...,tk on the axis in the orthogonal direction v. 
The dotted lines through these points represent measurement lines. A pair of points (in Fig. [TJ one red and 
one blue, shown in purple if they coincide) is placed randomly on each of the 4fc measurement lines. Since 
the geometry of the measurement lines is known, the position of each point can be described by a single real 
variable giving the location of the point on the measurement line. Therefore the position of all of the points 
can be described by a single vector variable z with 8fc real components. 

An initial guess for K is obtained by forming the convex hull of all 8k points, except those for which a 
pair coincides, i.e., the purple points. This is the convex polygon labeled P[z] in Fig. [T] The convex hull is 
computed using a standard algorithm as a subroutine. The reason for ignoring the purple points in taking 
the convex hull is that if a measurement line does not meet K , there must be some mechanism to eliminate 
the pair of points that lies on that line. In practice, a threshold is set so that a pair of points is eliminated 
if they become too close in the iterative optimization procedure to be described next. 

In order to improve the initial guess, the positions of the pairs of points on the measurement lines must be 
adjusted. This is effected by computing the sum, over all measurement lines, of the squares of the differences 
between the measured projection value for K and the corresponding projection value of P[z\. This least 
squares sum is the objective function in an optimization problem with 8k real variables and an optimization 
routine is used to drive the value of the objective function down to a minimum. The output of the algorithm 
is the convex polygon P^ = P[z] corresponding to the optimal vector z of these real variables. 

In [in] it is shown that for any finite set of directions for which the corresponding exact projections 
determine a convex object uniquely, the output Pk converges to K as k — > oo, even when the measurements 
are affected by Gaussian noise of fixed variance. Moreover, this remains true even if the optimization problem 
is not solved exactly, but only within an error Ek > 0, provided Ek — > as k — ¥ oo; see [HO p. 337]. In practice, 
the optimization problem involved is heavily non-linear. The f mincon function from Matlab's Optimization 
Toolbox was used, along with simulated annealing to improve performance. 




Figure 1: Illustration showing the basic principles behind the GKXR algorithm. 
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2.2.2 Unfiltered backprojection (U-FBP) 

U-FBP, as described here, can be viewed as a geometric method. The idea is to backproject the object's 
shadows (yielding a strip for each shadow) and to return the intersection of all of these strips. The returned 
object is then necessarily a convex polygon. This, in general, cannot be guaranteed for other common U-FBP 
variants that backproject projections instead of shadows. 

2.2.3 Modified Prince-Willsky (MPW) 

The modified Prince- Willsky algorithm [20] , a modification of the algorithm in [47] , reconstructs a convex 
object K from its support function hx. The function hx takes a direction (unit vector u) as input and 
returns a number that corresponds to the extent of K in direction u. To be precise, 

^if(u) = maxw T x, 

xEK 

where u T x denotes the inner product of u and x. 




Figure 2: Left to right: A convex set K in the plane; a set L which is not convex; the support function of 
K in one direction u; support function values in many directions describing the convex set K . 

Fig. [2] also indicates that K is completely determined by its support function values in all directions (this 
can be shown mathematically; see jTTJ Section 0.6]). Good approximations can already be obtained using a 
finite number of (suitably chosen) directions. 

The support function values of K in the two directions perpendicular to the projection direction can be 
easily determined from the data, because they correspond to the minimal (respectively, maximal) coordinates 
of the pixels in the data that record non-zero intensities. As data are available for many tilt angles, support 
function measurements are collected for different u vectors. These serve as input to the algorithm. 

So far, the algorithm is very similar to U-FBP. However, U-FBP tries to find an object that fits the 
noisy measurements perfectly. As with most inverse problems, this is usually not the best strategy. Fig. [3] 
illustrates the fact that noise may lead to inconsistencies in the data. 




Figure 3: Left: A convex set K and three exact measurements of its support function. Middle: An in- 
correct support function measurement leading to inconsistent data. Right: An incorrect support function 
measurement cutting away a big part of K. 

The MPW algorithm is designed to deal with noise. More precisely, for (Gaussian) noise affected mea- 
surements hi, . . . , h n of the support function of K in a finite number of directions ui, . . . , u n , the MPW 
algorithm solves a (linearly) constrained least-squares problem to obtain values yi, . . . ,y n , which are the 
support function values of the best-approximating set K*. 

With mild restrictions, the output of the algorithm converges as the number of shadows, affected by 
Gaussian noise of fixed variance, approaches infinity [21]. The implementation of MPW is somewhat more 
demanding than U-FBP, because a subroutine for solving quadratic programs is required (we use the Xprcss 
solver). The set K* , in our implementation, is obtained as an intersection of halfspaccs K* = f}™ = i{x '■ 
u Jx < Hi} via U-FBP. Additional MPW variants are discussed in [20] . 

2.2.4 2n-GON 

Again, shadows are taken as input data. We aim to reconstruct an object K that is known to be close 
to a regular 2n-gon, where n > 3 is known in advance. (In our experimental application, n = 3.) If the 
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assumptions are not fulfilled then the algorithm should exit without reconstruction. 

The length of the shadow of K for a given tilt angle 9 is commonly referred to as the width of K 
orthogonal to 9. Of course, this quantity can easily be computed from the input data. It is easy to see that 
if K is a regular 2n-gon, then the width of K, as a function of 9 € [0°, 180°), has exactly n local minima, 
corresponding to the tilt angles that project K along an edge direction of K . These n minima are thus 
(180/n)° apart and if they can be determined, then unfiltered backprojection (as implemented in U-FBP) 
from the corresponding n directions yields the regular 2n-gon K . Note that the shadows from two such edge 
directions of a regular 2n-gon K determine the minimum width and center of K and hence K itself. Also 
note that two such shadows are typically available from the data, because standard electron microscopes 
allow [0°, 120°) tilt ranges. 

The procedure also works when K is only close to a regular 2n-gon, as follows. Again, the idea is to apply 
unfiltered backprojection for the directions orthogonal to the tilt angles determined by some n minima of 
the width function of K that are only approximately (180/n)° apart. As before, some of the data for edge 
projecting directions might not be available if the data are not acquired over the full [0°, 180°) tilt range. 
If this is the case, then one needs to impose some assumptions on the missing data; see Fig. Eta), in which 
the hexagon and the parallelogram are indistinguishable from data in the (very limited) [0°,70°) tilt range 
shown. Another limitation of the 2n-GON approach is that we need to assume that the noise level in the 
(shadow) data is sufficiently low to allow determination of the minima of the width function. 

Here is a precise description of our implementation. We assume that data are available over a [0°,uj°) 
tilt range, where u> is fixed (typically to 140). We initially determine the local minima in [0°,u; ) of a 
polynomial curve of degree at least k > 2n that best fits (in the least-squares sense) the measured widths. We 
found in simulations that values of k around 2n + 5 give good and stable results. Fig. E^b) shows measured 
widths of a hexagon together with its best fitting polynomial curve. 




Angle (deg) 



Figure 4: (a) A regular hexagon (lightly shaded) and the parallelogram determined by width data from 
angles indicated by the arrows; (b) measured width as obtained from the HAADF STEM data (red) and 
best fitting (in the least-squares sense) polynomial curve of degree 11 (black). 

Suppose there are m > 2 minima ti < ti < ■ ■ • < t m , in degrees, within the [0°, u>°) range (otherwise, we 
stop without reconstruction). Unfiltered backprojection from S — {t±, . . . , t m } yields a (possibly degenerate) 
2m-gon, which, for ui = 180, approximates K. If u < 180, we proceed as follows. Let 

R= {i e {l,...,m- 1} : 180/n - 10 < \t i+1 - U\ < 180/n + 10} . 

If R = 0, we stop without reconstruction. If \R\ = n — 1, we reconstruct from the angles in R. Otherwise, 
1 < \R\ = r < n — 1 and we define a = max{i?}, b = a + 1, and T = {t[, . . . , t' n _ r }, where i'- = tb + 180i/n, 
i = 1, . . . , n — r. If there is an angle in T that is not in the [0°, 180°) range, we exit without reconstruction. 
Otherwise, we reconstruct from S U T. Here if t[ is outside the tilt range, we use unfiltered backprojection 
of the shadows for the angles t a and tb to obtain a parallelogram and set the shadow for t[ to be equal in 
length to that for tb, using the center of the parallelogram to position it correctly. 

We remark that this implementation reduces in the regular 2n-gon case to the method described at 
the beginning of this subsection. Note that 2n-GON can also be applied to data from a small number of 
tilt angles, since the best fitting polynomial curve may still approximate the local minima, although these 
minima might not be present in the data. We test this, among other things, in the next two sections. A 
theoretical analysis, however, remains outside the scope of this paper, as tolerable deviations from a regular 
2n-gon depend on several parameters such as the relative position of the missing wedge, the amount of noise 
in the data, the number of available projections, and the diameter of the 2n-gon. 

It is worth mentioning that as well as measuring widths orthogonal to 9, we can also measure widths 
parallel to 9, provided that projection data (and not only shadows) are available. We shall not discuss the 
performance of this variant here. 
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3 Experimental application 



We consider the task of reconstructing a nanowire from HAADF STEM data. Nanowires, small wires that 
are tens of nanometers in diameter and micrometers in length, are promising building blocks for future 
electronic and optical devices; see [351 142] . They are typically grown from a substrate and much research 
effort is being focused on understanding and controlling their growth mechanisms |13| . Electron tomography, 
as in various materials science applications, is rapidly developing into a powerful 3D imaging tool for studying 
these effects at the nanoscale 2TJ . 

With current technology, the tomographic data acquisition time for 140 projections of a single nanowire is 
about 2 hours when performed manually. This is currently a bottleneck preventing many in-situ experiments 
on short time scales and the imaging of multiple nanowires. Automated acquisition can lower this time, but 
a further reduction, which could be achieved if the required number of tilt angles can be reduced, is of 
paramount importance. Computation times for the reconstruction of a single 512 x 512 slice range from a 
few seconds to 3 minutes (for GKXR) on a standard PC. 

The particular nanowire in our experimental application is grown from pure InAs. Such nanowires are 
usually convex, or nearly so, in cross-sections perpendicular to the growth direction. In fact, most of these 
cross-sections for InAs nanowires are close to regular hexagons [52] . 

3.1 Experiment 

An HAADF STEM image series of an InAs nanowire was acquired using a probe-aberration-corrected FEI 
Titan 80-300 microscope operated at 300 kV with a probe convergence semi-angle of 18 mrad and an inner 
detector semi-angle of 100 mrad. The images were acquired over a total angular range of 139° with a 1° 
tilt increment. The original 2048 x 2048 pixel images, representing the projections for 2048 slices, were then 
binned by a factor of 4 to reduce noise in the reconstruction. The resolution after binning corresponds to a 
pixel size of 0.84 nm x 0.84 nm. 

Fig. [5ja-c) show a bright-field transmission electron microscopy (TEM) image and two HAADF STEM 
images, respectively, of the InAs nanowire specimen. An effect of non-linear projection intensities within a 
particular slice is visible in Fig. [5jd,e), as the measured projections (red graphs) do not exactly match the 
projections of the best-fitting shapes of the nanowire for this slice (black graphs). Twinning in the growth 
direction produces the "bee stripe" patterns visible in Fig. [5jb,c). 




Figure 5: Nanowire data: (a) bright-field TEM image of the InAs nanowire specimen used for tomography; 
(b,c) aligned binned HAADF STEM images taken from the tilt series of the wire at angles of 10° and 40°, 
respectively (slice 220 is indicated by the white line; the tilt axis is in the vertical direction through the 
center of the image); (d,e) measured projections of slice 220 (non-linear projection intensities) shown in red 
and ideal projections of slice 220 (linear projection intensities) shown in black, at angles of 10° and 40°, 
respectively. (The ideal projections were estimated from our reconstructions.) 
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3.2 Parameters for the algorithms 

The relaxation parameter A for SIRT and DART was set to 1. The SIRT algorithm was run for 50 iterations, 
based directly on the measured intensity values in the projection data. Afterwards, the reconstruction 
was thresholded at a value of 0.5, where represents the background and 1 represents the gray level of the 
nanowire. Each run of the DART algorithm consisted of 25 SIRT iterations to compute the starting solution, 
followed by 25 DART iterations, each of which included 10 SIRT iterations on the set of free pixels (i.e., 
the pixels near the boundary). The SIRT iterations within each DART iteration only take a fraction of the 
time of a complete SIRT iteration, as they are only applied to the free pixels. As for SIRT, DART employs 
a threshold at 0.5. A fixed fraction of 0.85 was used, meaning that 85% of all non-boundary pixels (selected 
randomly) are kept fixed in each DART step; the fixed fraction can be adapted to a specific noise level for 
more accurate results. We refer to [7] for details. 

The BART parameters, along with the code, are provided in Appendix A, Table [TJ 

For reconstructions with GKXR, the projection in each direction was measured at 40 equally spaced 
positions along an axis orthogonal to the direction and passing through the center of a circular window 
known to contain the object to be reconstructed (i.e., k = 40 in Fig.[l}. Simulated annealing was used, with 
a fixed cooling schedule. The projection data was pre-processed by a simple smoothing algorithm from digital 
signal processing, an averaging filter based on the IIR (infinite impulse response) design (see, for example, 
[31)1 Chapter 8]). More specifically, if the non-zero 512 projection measurements for a certain direction are 

. . . , y m , a smoothed set of measurements y[, . . . , y' m is produced by recursively defining 

y't = (yi + ?4_i)A 

for i = l,...,m, where we set tjq = 0. This smoothing procedure was iterated 50 times to obtain the 
smoothed data for input into the algorithm. 

The algorithms U-FBP, MPW, and 2?i-GON require shadows as input. Let Mg denote the 512 x 512 
matrix containing (row-wise) the projections of the nanowire slices for viewing angle 9. For this particular 
data set, we binarize the median filtered M$ (with a 1 x 3 window) using threshold T = 0, which is followed 
by a morphological opening [33l Chapter 15] with Matlab's structuring element strel (' diamond ' , 2). 
This morphological opening takes neighboring slices into account; the reconstruction of the object, however, 
proceeds slice by slice. 

3.3 Experimental results 

A full discussion of our results is only possible in connection with the simulations presented in the next 
section, particularly because there are no data for this nanowire obtained by an independent imaging method. 
However, we can make some general remarks. 

(i) It appears that the nanowire consists of two parts, the top part slightly narrower and rotated by about 
30° around the axis of the bottom part. The bottom part seems to consist of alternating twins, each 
around 40 nm thick, while one twin orientation dominates the top 80 nm. Except for the twinning (much 
less visible in GKXR and 2n-GON), these structures can be found in each of the seven reconstructions. 
However, the region around the frontal vertical edge of the reconstructions lies in the missing wedge 
and so might contain artefacts. Perhaps the most reliable indication of the presence of twinning is the 
fine periodic structure that appears along the left- and right-most vertical edges of the reconstructions. 
This periodic structure is too fine to be resolved by our current implementation of GKXR. 

(ii) BART and DART, followed by SIRT, appear to give the best (most accurate and precise) results. Note, 
however, that the software (Amira) that renders the 3D images in Fig. [6] smooths out isolated pixels 
due to the surface mesh that is built into the isosurface operation. Therefore these images might look 
smoother than those produced by these three algorithms without this post-processing. This is not the 
case for the other four algorithms, because they do not return objects containing isolated pixels. An 
explanation for the fuzzy boundary returned by GKXR is given in Section 14.31 

(iii) It seems possible to infer useful geometric parameters about the facet structure (diameters, angles, 
etc.) from GKXR, U-FBP, MPW, and 2n-GON, which use fewer data (shadows or fewer projections, 
as appropriate). See also Section |4*T41 

(iv) Algorithms U-FBP, MPW, and 2n-GON need segmentation of the shadows, while SIRT, BART, and 
DART need segmentation of the reconstruction. Depending on the reliability of the measured intensi- 
ties, one type of segmentation (possibly involving filtering and thresholding) might be favorable over 
another. 
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Figure 6: Reconstruction of the nanowire using different algorithms. Top-to-bottom and frontal views are 
shown in the first and second row, respectively. GKXR requires only 4 projections; U-FBP, MPW, and 
2n-GON reconstruct from shadows. A few of the top slices of the BART reconstruction are missing, because 
we put no effort into recovering the correct factor that relates the line integrals and the actual measured 
intensities. Slices are missing in the 2n-GON reconstruction if the algorithm detects that there is no hexagon 
in the corresponding slice. Nanowire orientations and their top-to-bottom views might vary as the 3D viewing 
points for the individual images have been manually selected. The 30 degree rotation of the end of the wire 
relative to the main part is more readily discernible in the frontal views (lower row) than in the top-to-bottom 
views (upper row). Images have been rendered using the Amira software with the isosurface operation and 
the compactif y option, but no explicit smoothing was applied. 

(v) It should be possible to improve the performance of GKXR by replacing the pre-processing via IIR 
with a least-squares fit to the projection data of a piecewise linear concave curve. Further improvement 
of the reconstruction quality of each of the presented algorithms might be achieved by incorporating 
the fact that the slices of the object are not independent from each other. The reconstruction from 
GKXR, for instance, might benefit from post-processing, such as a simple averaging process to smooth 
out the differences between consecutive slices. 

4 Computer simulations 

We tested the algorithms with four phantoms (i.e., simulated objects) under varying magnitudes of noise. 
As phantoms we used two regular hexagons, a slightly irregular hexagon, and a regular octagon. They are 
shown in Fig. [7] and henceforth are referred to as Phantoms 1-4. 

To compare phantoms and reconstructions we need to quantify the degree to which two shapes differ. 
(This is a central problem in computer vision.) Here we use two metrics that are frequently encountered in 
the literature, the symmetric difference metric and the Hausdorff metric. The symmetric difference distance 
between finite sets A and B of pixels (in our case, subsets of {1, . . . , 512} x {1, . . . , 512}) is defined by 

S S (A,B) = \(A\B)U(B\A)\, 

i.e., Ss(A, B) counts the number of mismatching pixels. The corresponding Hausdorff distance is defined by 

5 H {A,B) = max(6(A,B),5(B,A)) , 

where 

8(A,B) = max min I la — b\ |oo = max min max (\a\ — b\\, \a2 — M) 

a£A bEB {ai,a 2 )eA (6i,6 2 )eB 

and | • | denotes the usual absolute value of a real number. (We chose the norm for computational 
convenience; other choices are possible.) In other words, for S(A, B) we identify the pixel a £ A that is 
farthest from any pixel in B and return the distance from a to its nearest neighbor in B. Taking the 
maximum of 6{A, B) and 5{B, A) makes 8h symmetric in its arguments. Thus, roughly speaking, Sh (A, B) 
measures the extent to which each pixel of A lies near some pixel of B and vice versa. 
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Figure 7: Phantoms and reconstructions I. Column 1: Phantoms 1-4 (from top to bottom). Subsequent 
columns show difference images between a phantom and a typical reconstruction with 5140,10 at 50- noise 
obtained by SIRT (Column 2), BART (Column 3), DART (Column 4), GKXR (Column 5), U-FBP (Col- 
umn 6), MPW (Column 7), and 2n-GON (Column 8). Color scheme for difference images: White pixels 
belong to the phantom (P) and the reconstruction (R), red pixels to R \ P, and blue pixels to P \ R. These 
are 192 x 192 pixel images in which the 160 pixel boundary has been cropped. 

Now, let P denote the set of pixels of a phantom and R the set of pixels of a reconstruction. We refer 
to 6s(P,R) and Sh(P,R), or 5s and 8h for short, as reconstruction errors (measured in the symmetric 
difference metric and Hausdorff metric, respectively). 

The Hausdorff and symmetric difference metrics have different and somewhat complementary charac- 
ters. While the Hausdorff metric is sensitive to single pixel outliers but robust to boundary perturbations, 
the symmetric difference metric is robust to single pixel outliers but sensitive to boundary perturbations 
(cf. Fig. [71 particularly the DART reconstructions). However, as for any metric, two sets are equal if and 
only if they have zero distance. 

4.1 Data generation 

To avoid the "inverse crime" [TTJ Section 5.3] of using the same model for generating the data for testing the 
algorithms as the one used in their design, we generate the projections from higher resolution versions of the 
phantoms (cf. [31] [1]). More specifically, from 2048 x 2048 pixel versions of the phantoms we generate the 
projections and bin them by a factor of 4. We thus aim at reconstructing the 512 x 512 pixel versions of the 
phantoms as shown in Fig. [7J No further scaling of the projection values is introduced, i.e., the projection 
values give the number of pixels of the phantom that lie on the corresponding line. (Our implementation 
employs the Matlab command imrotate using the bilinear and crop option.) Comparing our results to 
those obtained from projection data generated by SNARK09 [12] . we could not find significant differences. 
This could be expected, because our projections (before binning) are generated from suitably high resolution 
phantoms. The noise model that is described next is applied to the projections generated from the high 
resolution phantoms. 

Taking a simplistic approach, we simulate Gaussian noise. Hence, we specify noise by one parameter 
a. We draw independent and identically distributed zero-mean, cr 2 -variance normal random variables p{i,j) 
for each projection angle i and projection pixel j that has non-zero intensity. The p(i,j)'s are added to 
the intensities of the corresponding projection pixels; negative intensities are set to zero. This approach 
simulates additive Gaussian noise on the non-zero intensities. The underlying assumption in our noise model 
of adding noise only to pixels with non-negative intensities is that statistical noise effects affect the recorded 
signal but still allow detection of the object's shadows. At least for our HAADF STEM data of the nanowirc, 
this seems to be a plausible model. In our simulations we consider the 0- noise (i.e., a = 0, noise- free) and 
50-noise {a = 50) cases. Here 50-noise seems to be in agreement with the intensity variations present in the 
experimental data in Section [3] 
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4.2 Parameters for the algorithms 

For our algorithms, the parameters used for all simulations were the same as for the nanowire slice recon- 
struction described in Section HOI The pre-processing of the data was simplified in the sense that for GKXR 
the pre-processing procedure to smooth the data was not used in the 0-noise case and for U-FBP, MPW, 
and 2n-GON the input shadows were obtained directly by thresholding the projections with T = 0. 

4.3 Simulation results 

We employed the following four sets of tilt angles: 5 18 o,x = {1°, 2°, 3°, . . . , 180°}, 5x40,1 = {1°, 2°, 3°, . . . , 140°}, 
5180,10 = {1°, 11°, 21°, . . . , 171°}, and 5i4o,xo = {1°, 11°, 21°, . . . , 131°}. Here a tilt angle corresponds to 
a clockwise rotation of angle around the vertical (0, Indirection. While 5x4o,i corresponds to our experi- 
mental setting, we included the other sets of tilt angles to compare the performances of the algorithms with 
respect to the missing wedge and the total number of available tilt angles. For the main simulations reported 
in Figs. [8] and [SJ reconstructions with GKXR were always performed using the four angles {1°, 28°, 91°, 118°}, 
chosen to be nearest to a theoretically ideal set within the range 1°-140°, namely angles corresponding to 
the vectors (0,1), (1,2), (1,0), and (2,-1) (cf. Section [2XT]). 

For each algorithm, each set of tilt angles, and each noise level we performed 100 reconstructions. The 
mean reconstruction errors 6s and Sh together with their standard deviations are shown in Fig. [5] and [51 
respectively. 

Wc hrst discuss typical difference images of reconstructions of each phantom, shown in Columns 2-8 of 
Fig. [H obtained from SIRT, BART, DART, GKXR, U-FBP, MPW, and 2n-GON with tilt angles 5 14 o,iO at 
50-noise. Note that this means reconstruction from only 14 projections, ten times fewer than is typically 
employed in ET (the 5x4o,i case). White pixels in the color scheme of our difference images correspond 
to correctly reconstructed pixels, red pixels belong to the reconstruction but not the phantom, and blue 
pixels belong to the phantom but have not been reconstructed. The reconstruction errors (6s] Sh) for 
Phantoms 1-4, respectively, are: 



SIRT: 


(1,482; 13) 


(1,562; 15) 


(1,367; 13) 


(1,594; 16) 


BART: 


(523; 3) 


(513; 3) 


(491; 4) 


(494; 4); 


DART: 


(873; 16) 


(697; 14) 


(843; 7) 


(955; 8); 


GKXR: 


(693; 4) 


(1,174; 7) 


(933; 6) 


(670; 5); 


U-FBP: 


(1,285; 12) 


(922; 10) 


(1,258; 11) 


(1,022; 10) 


MPW: 


(1,282; 12) 


(905; 10) 


(1,227; 11) 


(948; 10); 


2n-GON: 


(215; 2) 


(386; 3) 


(653; 4) 


(146; 1). 



The difference images illustrate several characteristics of the different algorithms. For instance, the algo- 
rithms that are specifically designed to reconstruct from only a few directions, GKXR and 2n-GON, yield 
small errors, both in the symmetric difference metric and the Hausdorff metric. (Recall that GKXR uses 
only four projections.) On the other hand, the effect of the missing wedge can be seen in the 2n-GON and 
MPW reconstructions, as these algorithms reconstruct from shadows. The 2n-GON algorithm cannot fully 
compensate for the missing wedge, because Phantom 3 deviates from a regular structure. The fuzziness of 
the object boundary in the BART and, particularly, the DART reconstruction is a common phenomenon for 
pixel-based reconstruction methods that do not employ filters. 

We now turn to a discussion of Figs.l8l fT0l We refer to 0-noise and 50-noise as moderate noise levels, while 
rr-noise with x > 200 is a high noise level. At high noise levels, the reconstruction quality, except perhaps 
for MPW and 2n-GON, becomes very poor. Therefore Figs. [5] and O only display results for moderate noise 
levels. Typical results for high noise levels are presented in Fig. 1101 and discussed in (|viii[) below. 

The simulations for 5x4o,x at 50-noise resemble the experimental conditions of the nanowire reconstruction 
presented in the previous section; the corresponding simulation results arc depicted by the yellow bars in 
Figs.|8jb) and[9jb). Experimental conditions with a reduced number of tilt angles would be represented by 
5x4o,io at 50-noise. The key inferences from the simulations can be summarized as follows. 

(i) For available projections over the whole 180° tilt range (as in the 5xso.i and 5xso,io-case) and moderate 
noise levels, we find that all algorithms yield good (accurate and precise) results (usually with Sh < 5 
and S s < 1000). 

(ii) BART, DART, and 2n-GON give the best results in the 5x4o,x case at moderate noise levels (the 
50-noise case resembles the experimental setup). 

(iii) 2n-GON performs well, even with missing wedge and few available projection directions. However, 
the object needs to be close to a regular 2n-gon, otherwise caution is required (see the results for 
Phantom 3). 
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Figure 8: Simulation results presented as mean reconstruction errors measured in the symmetric difference 
metric: (a) 5180,1: (b) 5*140,1, (c) Siso, 10, (d) Si4o,H> Bar colors indicate noise level, green for 0- and yellow 
for 50-noise. Black error bars represent standard deviation. GKXR requires only 4 projections; U-FBP, 
MPW, and 2n-GON reconstruct from shadows. 



(iv) In many cases, BART and DART give the best results. Particularly with DART, the quality of the 
reconstruction deteriorates significantly if fewer projections are available and if the noise level increases. 

(v) According to the simulations, the GKXR algorithm is among the better-performing algorithms in the 
missing wedge case at moderate noise (see Figs. [S^d) and^d)) and its outstanding feature is that it 
requires projection data only from four directions. To test the effect of changing the four directions 
used, we performed the same simulations with angles {21°, 41°, 61°, 81°}, all contained in a narrow 60° 
range. The results were worse, but not dramatically so. For example, for Phantom 1 at 50-noise, the 
mean (<5 5 ; S H ) errors rose to (1, 227.83; 8.48), compared to (925.84; 6.94) for angles {1°, 28°, 91°, 118°}. 
However, the nanowire reconstruction in Fig. [5] seems the worst produced by the algorithms used. To 
some extent, this is due to the fuzzy nature of the boundary depicted there. In fact, while the mean 
and standard deviations of the errors reported for GKXR in Figs. Hfb) and H^b) are fairly small, the 
distance between neighboring slices of the nanowire reconstruction can be significant. This is caused by 
the inherent stochastic nature of the algorithm, which uses simulated annealing. Possible improvements 
have been discussed in Section I3~2"1 

(vi) The problems for MPW and U-FBP are caused by the missing wedge, because data are simply missing 
from the shadows (see Fig. BJa)). If essential object features do not lie in the missing wedge, then 
good reconstruction results can be obtained even at high noise levels. In fact, MPW is among the 
better-performing algorithms if data are acquired over a 180° angular range. The U-FBP algorithm 
never outperforms MPW by a significant margin and in most cases the results for MPW are clearly 
superior to those reported for U-FBP. Inspection of the reconstructions shows that U-FBP tends to cut 
off the corners of the object, while MPW returns the correct shape up to a systematic overestimation. 
See also f|viii[) . 

(vii) SIRT is known for its noise suppressing features. It performs much better in highly limited data 
scenarios than transformation based techniques such as filtered backprojection |321 Chapter 7]. These 
findings are confirmed in our simulations. To avoid overloading the figures, however, we chose not to 
present the corresponding poorer results for those scenarios that we obtained by filtered backprojection. 
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Figure 9: Simulation results presented as mean reconstruction errors measured in the Hausdorff metric: (a) 
<Si80,i, (b) Si 4o,i, (c) Si8o,io, (d) Si4o,io- Bar colors indicate noise level, green for 0- and yellow for 50-noise. 
Black error bars represent standard deviation. GKXR requires only 4 projections; U-FBP, MPW, and 2n- 
GON reconstruct from shadows. (For clarity, we do not show the bars of the SIRT and DART results for 
Si4o,io with 50-noise if they extend beyond the value 15; their actual values range between 20 and 160.) 



While SIRT gives rather good results as measured by the Hausdorff metric for 180° angular range 
(Si8o,i, Si8o,io) at moderate noise levels, we observe from Figs. [5] and that the algorithm is typically 
outperformed by the other algorithms in the missing wedge case (Si4o,i, Si4o,io)- 

(viii) The reconstruction quality of all of the algorithms deteriorates at high noise levels. Fig. [TO] indicates 
that U-FBP, MPW, and 2n-GON seem to be less affected. (We show only results for Phantom 2 and 
Si 40,1 ; the other cases are similar.) This finding is somewhat expected, since these algorithms work 
with shadow data. 

(ix) There is a notable discrepancy between the performance of 2n-GON and GKXR in the Si 40.1 simulations 
with 50-noise and their performance, much worse in the case of GKXR, in the nanowire reconstructions 
shown in Fig. [6] Possible reasons for this discrepancy were given in ([m]) and (jvj). A detailed assessment, 
however, remains for future research. 

4.4 Determination of faceting 

A typical experimental aim when imaging geometric objects is to determine their faceting including minor 
facets and surface protrusion and roughness. Fig. 1111 shows two more complex phantoms along with recon- 
structions obtained by SIRT, BART, DART, GKXR, U-FBP, MPW, and 2?i-GON, respectively. We refer 
to these as Phantoms 5 and 6; reconstructions for demonstration purposes have been performed for Si4o,io 
at 50-noise (the data generation was described in Section l4~Tj) . 

The reconstruction errors (5s; Sh) for Phantoms 5 and 6 are, respectively: For SIRT, (1,501; 13) and 
(1,742; 12); for BART, (460; 3) and (692; 9); for DART, (842; 11) and (769; 14); for GKXR, (934; 5) and 
(1,744; 11); for U-FBP, (852; 10) and (1,199; 8;); for MPW, (791; 10) and (1,243; 8); and for 2n-GON, 
(379; 5), and (1,443; 12). 

The results in Fig. [TT] can be summarized as follows. 

• SIRT: Minor facets of Phantom 6 are not reconstructed; the boundaries are very fuzzy; additional 
imaging tools for edge extraction are required to evaluate facet angles; the inner angles of the major 
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Figure 10: Simulation results for Phantom 2 with 5140,1 including 200-noise, presented as mean reconstruc- 
tion errors measured in (a) the symmetric difference metric and (b) the Hausdorff metric. Bar colors indicate 
noise level, green for 0-, yellow for 50-, and red for 200-noise. Black error bars represent standard deviation. 
GKXR requires only 4 projections; U-FBP, MPW, and 2n-GON reconstruct from shadows. 
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Figure 11: Phantoms and reconstructions II. Column 1: Phantoms 5 (top) and 6 (bottom). Difference images 
between a phantom and a typical reconstruction with Si4o,io at 50-noise obtained by SIRT (Column 2), 
BART (Column 3), DART (Column 4), GKXR (Column 5), U-FBP (Column 6), MPW (Column 7), and 
2n-GON (Column 8). Color scheme for difference images: White pixels belong to the phantom (P) and the 
reconstruction (R), red pixels to R \ P, and blue pixels to P\R. These are 192 x 192 pixel images in which 
the 160 pixel boundary has been cropped. 



facets are reconstructed with an accuracy of about 12 degrees. 

• BART: All minor facets are reconstructed; artefacts in the missing wedge region might be misclassified 
as minor facets; additional imaging tools for edge extraction arc required to evaluate facet angles; the 
inner angles of the major facets are reconstructed with an accuracy of about 2 degrees. 

• DART: Several minor facets are reconstructed; a few minor facets in the missing wedge region are 
not resolved; boundaries can be fuzzy and additional imaging tools for edge extraction are required to 
evaluate facet angles; assuming that only isolated pixels are removed, we find that the inner angles of 
the major facets are reconstructed with an accuracy of about 6 degrees. 

• GKXR: Few minor facets are reconstructed; the inner angles of the major facets are reconstructed with 
an accuracy of about 8 degrees (cf. Fig. [7]). 

• U-FBP and MPW: Minor facets of Phantom 6 are not reconstructed; the inner angles of the major 
facets outside the missing wedge are reconstructed with an accuracy of about 2 degrees; artificial major 
facets appear in the missing wedge of Phantom 5. 

• 2n-GON: Minor facets of Phantom 6 and two minor facets of Phantom 5, are not reconstructed; the 
inner angles of the major facets are reconstructed with an accuracy of about 1.5 degrees. 

These results demonstrate the potential of the geometric approach. A detailed study, however, must be left 
for future research, since the reconstruction quality depends on several parameters, such as the particular 
type of noise, the misalignment of the projections, and the relative position of the object and the missing 
wedge. 

As a final comment, we remark that a large I/2-norm of the residual Ax — b can be seen as an indicator 
that the reconstruction cannot be trusted (notation as in Section 12. 1 . 1[) . For the GKXR, MPW, and 2n-GON 
reconstructions in Fig. [TTJ for instance, we find that \\Ax — b\\ is at least three times larger for Phantom 6 
than for Phantom 5. 
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PICTURE TEST 

PROJECTION REAL 

MODE LOWER = 0.0 

MODE UPPER = 1.0 

STOP ITERATION 5 

EXECUTE AVERAGE ART CONTOUR 

BART on hexagon image 

0.5 1 

1 

ART3 relaxation constant 0.3 
CONSTRAINT BART 



STOP ITERATION 1 
BASIS PIXEL 

EXECUTE CONTINUE ALP1 SMOOTH 
Smoothing of BART recon 
2 111 
1 

BASIS BLOBS 

EXECUTE CONTINUE ALB1 CONTOUR 
Contouring of smooth BART recon 
0.5 1 
1 



END 



Tabic 1: SNARK09 code (input file) for the BART algorithm. Left: main routine. Right: filtering routine. 
Seven copies of the filtering routine should be placed between the "*" and "END" of the main routine. 

5 Conclusion 

We have applied five algorithms from the mathematical fields of geometric and discrete tomography to 
reconstruct homogeneous objects from electron tomography data. Our results demonstrate that the choice of 
reconstruction algorithm should be based on the specific reconstruction task at hand. None of the algorithms 
considered is better for all reconstruction tasks than the others. The main features of the algorithms that we 
introduced here in the ET context can be summarized as follows: BART and DART reconstruct homogeneous 
objects from projections; GKXR reconstructs convex objects from only four projections; MPW reconstructs 
convex objects from shadows and compensates for noise; and 2n-GON reconstructs objects that are close to 
regular 2n-gons from (possibly few) noisy shadows. 
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